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SUMMARY 

The use of various pressure correction algorithms originally developed 
for fully staggered meshes can yield a velocity-pressure decoupled solution 
for collocated meshes. The mechanism that causes velocity-pressure 
decoupling is identified in this paper. It is shown that the use of a 
partial differential equation for the incremental pressure eliminates such 
a mechanism and yields a velocity-pressure coupled solution. Example flows 
considered are a three-dimensional lid-driven cavity flow and a laminar 
flow through a 90 degree bend square duct. Numerical results obtained using 
the collocated mesh are in good agreement with the measured data and other 
numerical results. 
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Nomenclature 


*>h 

(i ,m,n) 

n j 

P 

Re 

U j 
X j 

V 

Superscripts 

n 

nb 


coefficients for incremental velocity, (i-=u, v, or w) 
hydraulic diameter of the curved duct 
index for mesh 

outward normal vector, nj-tn^ny ,n z } 

pressure 

Reynolds number 

velocity component, Uj«(u,v,w} 

cartesian coordinates Xj-{x,y,z} 

kinematic viscosity of fluid 

curvilinear coordinates £j-(£>r;,0 


iteration count 

neighboring grid points; (i+l,m,n), (i-l,m,n), (i,m-l,n), 
(i,m+l,n), (i,m,n-l) and (i.m,n+l) 

current value 

incremental (or corrective) value 


Subscripts 

i,j index for spatial coordinates, i— { 1,2,3} and j={l,2,3} 

Mathematical symbol 

Y summation 
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INTRODUCTION 


Finite volume methods based on various pressure correction algorithms 
are still the most widely used numerical methods for incompressible flows. 
The fully staggered mesh^*^ is used in most of these numerical methods to 
avoid velocity-pressure decoupling. However, application of the fully 
staggered mesh for flows with complex geometries is not straight forward, 
especially for three-dimensional flows. For example, the fully staggered 
mesh can not be used to solve flows inside a 90 -degree bend duct, see 
Maliska and Raithby^ for more details. A number of grid layouts have been 
proposed and are in use to solve the Navier- Stokes equations defined on 
complex geometries. In Maliska and Raithby^, the standard fully staggered 
mesh is extended in such a way that the three velocity components are 
located at all grid points except at the pressure grid point. In this case, 
the number of degrees of freedom for velocity is tripled while the degrees 
of freedom for pressure remains the same as that of the original fully 

staggered mesh. In Vanka et al.^ and Chen^ , the velocities are located at 

f 

the same grid points and the pressure is located at the centroid of the 
cell formed by the four adjacent velocity grid points (for the 
two-dimensional case). Vanka et al.^ mentioned that it was not easy to 
obtain convergent solutions due to velocity-pressure decoupling. This grid 
layout is denoted as the pressure -staggered mesh (ps-mesh) in the following 
discussion. In Chen^ , velocity-pressure decoupling Is eliminated by using a 
non- conforming domain for the mass imbalance calculation. A collocated mesh 
is used in Rhie^, Dwyer and Ibrani^, Peric et al.® and Majumdar^. In Rhie^, 
velocity-pressure decoupling is prevented by including a fourth order 
artificial dissipation in the pressure correction equation, while in Dwyer 
and Ibrani^ , the same purpose is achieved by evaluating the incremental 
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velocities at mid-sides of the control volume. In Peric et al.^ and 
Majumdar 9 , velocity-pressure coupling is achieved by interpolating the 
pressure gradient term differently from the rest of the discrete momentum 
equation. This interpolation scheme is called "momentum interpolation” in 
Majumdar 9 . The momentum interpolation scheme becomes similar to the present 
pressure correction algorithm in the case of an orthogonal mesh, and is 
discussed in detail in the following sections. From the above discussion, 
the technique for coupling the velocity and pressure forms the basis for 
each numerical method. In other words, the grid layout that can be used is 
limited by the velocity-pressure coupling mechanism to be adopted. 

In the case of a collocated mesh, the velocities at odd-numbered grid 
points are driven by pressure at even-numbered grid points, and vice versa. 
Thus the momentum equation can not sense the presence of a checker -board 
type pressure field 1 . Thus it is believed that a collocated mesh yields a 
velocity-pressure decoupled solution. However, a number of numerical 
results obtained using a collocated mesh have been reported In recent 
years^ -1 ^. The successful calculations of various flows obtained using the 
collocated mesh suggest that velocity-pressure decoupling may be caused by 
various pressure correction algorithms rather than by the collocated mesh 
itself. A re-examination of the velocity-pressure decoupling mechanism in 
collocated mesh solutions and a partial differential equation for the 
incremental pressure that can be used for various grid layouts are 
presented In this paper. 


ANALYSIS 

The incompressible laminar flow equations are given as; 
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( 1 ) 


du-j 

— - 0 . 



where the subscripts i and j denote each coordinate direction, and the 
repeated indices imply summation over the index unless otherwise stated. 

In the pressure correction methods^-* the velocities and the pressure 
are decomposed as; 

Ui - + U' it (3) 

P - P* + P' (4) 

where the superscript * denotes the current velocities which may not 
satisfy the conservation of mass equation yet. The discrete momentum 
equation for u^-velocity can be written as; 

dp 

(2 ,m , n)u^ (i , m ,n) - £ ^i u i^ ~ + S^ v , no sum on i (5) 

nb 

where A^ is the coefficient of the u^-veloclty at velocity grid point 
(i,m,n), Si v is the load vector originating from the curvilinear grid 
structure, the superscript nb denotes the neighboring grid points, and the 
pressure gradient is left in continuous form deliberately. The major 
difference between the present pressure correction algorithm and the other 
algorithms originates from leaving the pressure gradient term in continuous 
form in deriving the pressure correction equation. In fact, 
velocity-pressure decoupling that occurs when various pressure correction 
algorithms are used for a collocated mesh or a pressure-staggered mesh is 
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caused by using a discretized pressure gradient in eq. (5) in deriving the 
pressure correction equation. The details are discussed later in this 
section. Substituting eqs. (3-4) into eq. (5) yields; 


9(P*+P') 


Ai( u i +u^') - £ {A^(u£ +U£ 1 ) ) — + Si v , no sum on I (6) 

nb 


where the grid point index (i,m,n) has been deleted in eq. (6) for 
convenience. The discrete -momentum equation based on the current flow 
variables which may not satisfy the conservation of mass equation can be 
written as; 


. . - r 

A-^u^ - Y, ( A i u i ) + $i V > no sum on ^ 

nb dx£ 


Subtracting eq. (7) from eq. (6), disregarding the summation over the 
neighboring grid points and the load vectors, yields; 


1 dp' 


, no sum on l 


A i $ x i 


Substituting eq. (3) into eq. (1) yields; 


Formally substituting eq . (8) into eq. (9) yields the elliptic partial 
differential equation for the incremental pressure; 


a f 1 ap' 


axj l a j axj 
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where the last terra in eq. (10) represents the mass imbalance. Applying the 
standard finite volume procedure to eq. (10) yields a system of equations 
for the incremental pressure. In this case, the value of Aj at the 
interface of the pressure control volumes is obtained by linearly 
Interpolating the adjacent nodal values. As all the central -differenced 
finite volume equations for self-adjoint second order elliptic partial 
differential equations are strongly diagonally dominant, the present 
discrete pressure correction equation is strongly diagonally dominant. 

The mechanism that leads to a velocity-pressure decoupled solution, 
when various pressure correction algorithms are used for collocated or 
pressure - staggered meshes, is described below. In various pressure 
correction algorithms, the discrete pressure correction equation is 
obtained by integrating eq. (9) over a pressure control volume using a 
discretized pressure gradient in eqs . (6-8). This procedure yields a 
27 -diagonal system of equations for the incremental pressure for a 
three-dimensional case. This discrete pressure correction equation is not 
diagonally dominant, and thus, the mass imbalance for a particular pressure 
grid point tends to correct only the pressure of the adjacent pressure grid 
points. Velocity-pressure decoupling occurs in such a case. In case a 
discretized pressure gradient is used in deriving the pressure correction 
equation, the discrete system of equations can be made diagonally dominant 
by using a fully staggered mesh or by using a special interpolation 
scheme® j n t ^ e momen tum interpolation scheme of Peric et al.® and 
Majumdar^, the pressure gradient in eq. (7) is expressed using the pressure 
at two grid points located across the surface of a control volume and the 
other terms are linearly Interpolated. The discrete pressure correction 
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equation obtained using the momentum interpolation scheme is also strongly 
diagonally dominant. For an orthogonal mesh, the present pressure 
correction algorithm and the momentum interpolation scheme® yield similar 
seven- diagonal systems of equations for incremental pressure. The pressure 
correction algorithm of Peric et al.® and Majumdar^ always yields a 
seven- diagonal system of equations for incremental pressure regardless of 
grid skewness, whereas, the present pressure correction algorithm yields a 
27 -diagonal system of equations for arbitrary skewed mesh. Through a number 
of numerical tests using the collocated mesh and the pressure -staggered 
mesh-^, it was found that retaining the off-diagonal terms (which originate 
from the grid skewness) in the discrete pressure correction equation 
accelerates the convergence rate slightly compared with disregarding the 
off-diagonal terms. However, the advantage of slightly faster convergence 
rate is overshadowed by the disadvantage of the increased memory size to 
store the off-diagonal terms, especially for three-dimensional flows. Thus 
the present method and that of Peric et al.® share similar numerical 
characteristics. In the momentum interpolation scheme, the entire residual 
obtained by subtracting eq . (7) from eq . (6) is retained in eq . (8.), and 
the residual (other than the incremental pressure gradient term) also drive 
the incremental pressure. Thus the numerical result may depend upon the 
under- relaxation parameter unless an improved momentum interpolation scheme 
is adopted^. On the other hand, the present method does not yield a 
numerical result that depends on the under -relaxation parameter as eq. (10) 
clearly indicates that the incremental pressure is driven only by the mass 
imbalance. Also the present pressure correction algorithm can be most 
easily extended to solve compressible flow by Including convective 
incremental pressure terms into eq . (10), see Ref. 11 for details. 
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The numerical procedures are briefly described below. The pressure 
gradient in the momentum equation, eq. (5), for 
three-dimensional curvilinear coordinates, is obtained as 
dp dp d $ j 

(i,m,n) ( 11 ) 

dxj_ 3£j 3x^ 

where , for example , 

dp 

— (i.n.n) - {p(i+l/2 ,m- 1/2 ,n- 1/2) + p(i+l/2 , m+1/2 ,n-l/2) 

+ p(i+l/2,m+l/2,n+l/2) + p(i+l/2 ,m-l/2 ,n+l/2) }/4A| 

- {p(i-l/2 ,m-l/2 ,n-l/2) + p(i-l/2 ,ra+l/2 ,n-l/2) 

+ p(i- 1/2 ,m+l/2 ,n+l/2) + p(i- 1/2 ,m-l/2 ,n+l/2) }/4A£ 

and p(i+l/2 ,m+l/2 ,n+l/2) is obtained by averaging the eight adjacent nodal 
pressures, see Fig. 1 for the grid index. The pressure gradient equation, 

eq. (11), is the same as that used in the pressure - staggered mesh^ and is 

also compatible with that proposed by Maliska and Raithby^. The same 
pressure gradient equation is used in evaluating the incremental velocities 
at the grid points, eq. (8), with the pressure in eq. (11) replaced by the 
incremental pressure. The power law upwind differencing scheme-^ is used for 
the momentum equation, and the discrete momentum equations and the discrete 
pressure correction equations are solved by the 
Tri- Diagonal -Matrix -Algorithm^ . 

NUMERICAL RESULTS 

Numerical results for a three-dimensional lid-driven cavity flow^ and 
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a laminar flow through a 90° bend square duct 1Zf , obtained using the 
collocated mesh, are presented below. 

The lid- driven cavity flow at Reynolds number (based on the transverse 
velocity and the length of the lid) 3200 is considered. The symmetric half 
of the flow domain is discretized by 78x67x67 grid points in x- , y- , and 
z-coordinate directions, respectively. The calculated transverse and 
vertical velocity profiles along the center lines of the symmetry plane are 
shown in Figs. 2- (a) and 2-(b), respectively. It can be seen in these 
figures that the present numerical results compare favorably with the 
measured data and the other numerical results obtained by Freitas et al. 1 ^ 1 
using the SIMPLE algorithm 1 and those obtained by Kim 10 using the 
pressure-staggered mesh 10 . 

A laminar flow through a 90° bend square duct 1 ^ Is considered below. 
The flow geometry is schematically shown in Fig. 3. The Reynolds number 

based on the hydraulic diameter of the duct (E^-0.04 m) and the bulk 

velocity is approximately 800. The upstream boundary is located at 10D h 

upstream of the curved section and the exit boundary is located at lOD^ 

downstream of the curved section. The symmetric half of the flow domain is 
discretized by a 92x18x33 grid points in x- , y- , and z-coordinate 
directions, respectively. The velocity profile of a fully developed square 
duct flow Is prescribed at the inlet boundary and the vanishing gradient 
boundary condition is used at the exit boundary. The calculated 
longitudinal velocity profiles are shown in Fig. 4. The numerical results 
obtained using the collocated mesh are almost the same as those obtained 
using the pressure- staggered mesh. It can be seen in the figure that the 
present numerical results are In good agreement with the measured data and 
the numerical results by Humphrey et al. 1 ^ The pressure distributions along 
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the inner and outer walls at the symmetric plane are shown in Fig. 5. The 
nodal pressure obtained using the collocated mesh exhibits mild wiggles. 
However, the averaged pressure at the cell interface varies smoothly as 
shown in Fig. 5. Again, the numerical results obtained using the collocated 
mesh are in good agreement with the other numerical results^ The 
projected velocity vectors on planes very close to the side and the bottom 
walls (0.0050^ away from each wall) are shown in Figs. 6-(a) and 6-(b), 
respectively. These planes are located so close to the wall that the 
velocity component normal to each plane is by far smaller than the velocity 
component shown in the figure. Due to the strong adverse pressure gradient, 
a small reversed flow region is formed near the corner wall and the 
reversed flow region extends up to 0=37° as can be seen in these figures. 
The numerical results obtained using the collocated mesh are essentially 
identical to those obtained using the pressure- staggered mesh, see Ref. 10 
for more details . 


DISCUSSION 

The use of various pressure correction algorithms^ » ^ for the 
collocated mesh yields ill-conditioned (e.g., very small diagonal entry) 
discrete pressure correction equation. In such a case, the mass imbalance 
at a particular pressure node corrects only the pressure of the adjacent 
pressure grid points, and a converged solution can not be obtained. On the 
other hand, applying the standard finite volume numerical procedure to the 
partial differential equation for incremental pressure yields a strongly 
diagonally dominant discrete pressure correction equation and a physically 
consistent numerical result is obtained. The present study and those by 
other investigators^'^-® indicate that velocity-pressure decoupling is 
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caused by an ill-conditioned discrete pressure correction equation. 

The present pressure correction algorithm and the "momentum 
interpolation" schemes of Peric et al . 8 and Majumdar 9 share similar 
numerical characteristics in the sense that these methods yield strongly 
diagonally dominant systems of discrete pressure correction equations. On 
the other hand, the present pressure correction algorithm does not involve 
a mixed interpolation of the discrete momentum equation as used in the 
momentum interpolation schemes 8 * 9 . Majumdar 9 pointed out that the use of a 
mixed interpolation for the discrete momentum equation may yield a 
numerical result that depends on the under- relaxation parameter and 
proposed an improved momentum interpolation scheme that yields a unique 
solution. The present pressure correction algorithm does not yield a 
numerical result that depends on the under -relaxation parameter, and it can 
also be used for pressure- staggered mesh^ and compressible flows^. 
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(a) On symmetric plane (y/D h = 0.5). 

(b) On quarter plane (y/D h =0.25). 

Figure 4. — Velocity profiles for curved duct flow. 
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Figure 5, — Pressure distributions on the wail. 



(a) (£, {) plane at yDh = 0.005 


<<<<<<<< <<<<<<o>>>>>> >>>>>>>>>>>>>>>>> 

v v vr t* tvw :>■?’>>>>>>>>>>>>>>>>>>>>>> 

"7 A 7777»>> >>>>>>>>>>>>>>>> 

A^A'1177777 77 7 7 7>7>7>>>>»>>»> 
AAAAA ill'f 1 tr 

^77 7^777^77?-7»>^»»>^> 
>>777711^77^777 T ^ ^ ^ ->* > > 

>>>7^7'7'7'f'T'7 7777 ^ ^ ^ ^ » 

>>>77777777 777> >>> >»>>> 

>>>>7777»>»>> > >-»>>»>> 


(b) (i ti) plane at £/Dh = 0.005. 


Figure 6. — Velocity vectors. 
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